Collapse to Black Holes in Brans-Dicke Theory: I. Horizon Boundary 
Conditions for Dynamical Spacetimes 

Mark A. Scheel 

Center for Radiophysics and Space Research 

and Department of Physics, 
Cornell University, Ithaca, New York 14853 

Stuart L. Shapiro and Saul A. Teukolsky 

Center for Radiophysics and Space Research 

and Departments of Astronomy and Physics, 

Cornell University, Ithaca, New York 14853 

November 1994 



ABSTRACT: We present a new numerical code that evolves a spherically symmetric configu- 
ration of collisionless matter in the Brans-Dicke theory of gravitation. In this theory the spacetime 
is dynamical even in spherical symmetry, where it can contain gravitational radiation. Our code is 
capable of accurately tracking collapse to a black hole in a dynamical spacetime arbitrarily far into 
the future, without encountering either coordinate pathologies or spacetime singularities. This is 
accomplished by truncating the spacetime at a spherical surface inside the apparent horizon, and 
subsequently solving the evolution and constraint equations only in the exterior region. We use 
our code to address a number of long-standing theoretical questions about collapse to black holes 
in Brans-Dicke theory. 



I. INTRODUCTION 



In recent years, there has been renewed interest in scalar-tensor theories of gravitation. One reason is 
that these theories are important for cosmological inflation models[l], in which the scalar field allows the 
inflationary epoch to end via bubble nucleation without the need for fine-tuning cosmological parameters 
(the "graceful exit" problem). In addition, scalar-tensor gravitation ("dilaton gravity") arises naturally from 
the low-energy limit of superstring theories[2, 3]. Finally, with the construction of LIGO, it may be possible 
to test scalar-tensor theories to high precision [4] by looking for monopole and dipole gravitational radiation 
from astrophysical sources. 

Quite apart from their potential physical significance, scalar-tensor theories play another very useful role: 
they provide an ideal laboratory for testing new algorithms for numerical relativity. In general relativity, 
numerical methods for treating spacetimes containing gravitational radiation require at least two spatial 
dimensions, since a time-varying quadrupole moment is needed to produce gravitational waves. In scalar- 
tensor theories, one can study many of the same strong-field phenomena that occur in general relativity, 
including gravitational radiation and dynamical black holes, while still working in spherical symmetry. 

We have developed a numerical code that solves the coupled matter and gravitational field equations 
for the evolution of a spherically symmetric configuration of noninteracting particles in Brans-Dicke[5] grav- 
itation, the simplest of the scalar-tensor theories. We use this code to study gravitational collapse to a 
black hole in Brans-Dicke theory. This process has been discussed extensively in the literature[6], but these 
studies have been limited to addressing the final state of the black hole after collapse, or have used linearized 
approximations of the field equations. Other than an early simulation by Matsuda and Nariai[7], it is only 
very recently[4] that this process has been calculated in any detail. 

In constructing numerical models of gravitational collapse in Brans-Dicke theory, we have been forced 
to address the same difficulty that has plagued the field of numerical relativity for the last 30 years: how 
does one handle the spacetime singularity at the origin that inevitably develops during the formation of a 
black hole? 

The traditional approach has been to utilize the "many-fingered time" gauge freedom of general relativity 
to avoid the singularity altogether. Specifically, one chooses coordinates such that the passage of proper time 
grinds to a halt near the origin before the singularity appears, while weak-field regions of spacetime farther 
from the origin evolve farther into the future. This singularity-avoiding (SA) method works well for short 
times, but eventually pathologies develop in the transition region between the "frozen" interior and the 
"evolving" exterior. These typically take the form of steep gradients or spikes in the metric functions, and 
will eventually cause the numerical code to crash [8]. Countermeasures such as increasing the grid resolution 
produce little improvement because the pathologies increase exponentially with time. 

Our solution to this problem is to use an apparent horizon boundary condition (AHBC) method after 
the formation of a black hole. The basic idea of this approach is to truncate a black hole spacetime at a 
surface inside the apparent horizon (AH) that cannot causally influence the exterior. One then discards the 
singular interior entirely, and only evolves the physically relevant exterior. 

Seidel and Suen[9] have implemented this idea in general relativity for the case of a black hole with a 
Klein-Gordon field in spherical symmetry. Their method involves a coordinate system that is locked to the 
AH, so that the coordinate speed of radially outgoing light rays inside the AH is negative. This enables them 
to use a causal difference scheme, similar to the Causal Reconnection Scheme of Alcubierre and Schutz[10], 
to solve evolution equations in such a way that information does not escape from the black hole, and no 
explicit boundary condition is needed on the AH. 

Our AHBC method is different from that of Seidel and Suen. Although we also use a coordinate system 
that is locked to the AH, we solve the wave equation for the Brans-Dicke scalar field using an implicit 
differencing scheme motivated by the work of Alcubierre[ll]. In addition, we solve for the metric variables 
using the elliptic constraint equations rather than the evolution equations. We obtain the required boundary 
conditions for this approach from asymptotic flatness, properties of the apparent horizon, and by solving a 
single evolution equation only on the AH, as explained in Section IV. 



With our code we are able to follow accurately Brans-Dicke collapse to black holes and the associated 
generation of monopole gravitational waves. We are able to integrate the equations to arbitrarily late times 
into the future, when all of the radiation has propagated out to large distances and the central black hole 
has settled down to final equilibrium. Using our code we are able to resolve a number of long-standing, 
theoretical questions about collapse in Brans-Dicke theory. We are also able to refine a promising technique 
for evolving black hole spacetimes with radiation by integrating only in the observable regions of spacetime. 

This paper is primarily concerned with numerical methods. The reader not interested in the numerical 
details should read section II and then skip to Paper II[12], in which we discuss Brans-Dicke gravitation in 
more detail, and we show how black holes formed in Brans-Dicke theory behave differently than those in 
general relativity. 



II. BASIC EQUATIONS 

A. Brans-Dicke Theory 
The action for Brans-Dicke gravitation is [5] 



where the Lagrangian density is 



c BD (-g) 1,2 d 4 x, (2.1) 



c BD = g ab R ab <f> + ^c-^ g ab d a <f>d b <f>. (2.2) 



The coupling constant to is dimensionless, and the scalar field <f> has dimensions of G _1 , where G is Newton's 
gravitational constant. The Lagrangian density C for matter and nongravitational fields depends on the 
metric g ab but not on <f>. The Ricci tensor R a i, is related to the metric in the usual way. In general relativity, 
the third term in Eq. (2.2) is absent, and one sets <f> = G _1 . In other scalar-tensor theories[13], the Lagrangian 
is more complicated. In particular, the coupling parameter to can be a function of <f>. 

In general relativity we are free to choose our units of mass and time such that G = c = 1. In Brans- 
Dicke theory, the inverse of the scalar field (f) -1 plays the role of G, so multiplying <f> by a global scaling 
factor changes the unit of mass. We choose units such that 

c=</> 00 = 1, (2.3) 

where (f)^ is the value of the scalar field far from any sources. 

Variation of Eq. (2.1) with respect to g ab and <f> yields the Brans-Dicke field equations, which can be 
written in the form 

□*=-^, (2.4) 

Y 3 + 2w' K J 

G ah = SttT^, (2.5) 



where 



87rT ab 



8irT ah + - (V^V^ - iffajV^V^) + V a V^ - g ai U<< 



(2.6) 



Here V denotes covariant differentiation with respect to the metric g a b, □ is the covariant Laplacian V a V a , 
and G a i, is the usual Einstein tensor. The symmetric tensor T a b is the energy-momentum tensor for matter 
and nongravitational fields, and T is its trace: 



T al = g al L-2— b , (2.7) 



T = T\=T ai g ah , (2.8) 

V b T ab = 0. (2.9) 

Although we have written the field equations (2.5) in a form that resembles Einstein's equations, we 
emphasize that Brans-Dicke gravitation is not the same as general relativity with a Klein-Gordon scalar 
field. The difference is the factor of <f> in the first term in the Lagrangian density (2.2), which leads to 
second derivatives of <f> in the field equations (2.5). Physically, this manifests itself as a violation of the weak 
equivalence principle for massive bodies[14]. This is discussed in more detail in Paper II. 

Notice that the matter stress-energy is conserved (Eq. 2.9), even though T ab is not equal to 87rG ai 
(the quantity VjT' ai also vanishes because of the Bianchi identity VjG ai = 0). As a result, the equations of 
motion of matter do not involve the scalar field; test particles move on geodesies of the metric. Notice also 
that it is the trace of the matter stress-energy tensor T a b, not the effective tensor T' ab , that appears in the 
wave equation (2.4). 

In vacuum, <f> = constant is a solution of Eq. (2.4) for any to. In this case, T' ah reduces to T a b and the 
field equations (2.5) become Einstein's equations. Therefore, any vacuum solution of Einstein's equations is 
also a vacuum solution of the Brans-Dicke equations. In addition, many (but not all) Brans-Dicke solutions 
with |w| — ► oo have <f> — ► constant and obey Einstein's equations. It is therefore said (but not rigorously 
correct — See Paper II) that Brans-Dicke theory reduces to general relativity in the limit |w| — ► oo. 

B. Equations (2.5) in (3+1) form 

Adopting the usual ADM[15] (3 + 1) decomposition, we introduce a set of spacelike hypersurfaces, or 
time slices, and a timelike vector field n a normal to these hypersurfaces. We adopt the convention that the 
4-metric g a b has signature ( \- ++), so that 

n a n a = -l. (2.10) 

Spatial distances on a particular time slice are measured by the 3-metric j a b, defined by 

jab = g a b + n a n h . (2-11) 

The extrinsic curvature K a b describes the rate of change of the three-metric along n a (the "time" direction): 

K ai = -\£ nlai . (2.12) 

Here £ denotes a Lie derivative. The field equations are split into spatial constraints that relate K a b and 
jab on each time slice, and first-order (in time) evolution equations that take K a b and j a b from one slice to 
the next. 

We work in spherical symmetry, and we choose the maximal time slicing and isotropic radial coordinate 
conditions. This gauge is defined by the isotropic line element 

ds 2 = -(a 2 - A 2 j3 2 )dt 2 + 2A 2 j3drdt + A 2 (dr 2 + r 2 dtt 2 ) (2.13) 

and the maximal slicing condition 

K = K )t = 0. (2.14) 
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Here K is the trace of the extrinsic curvature K a b- The 3-dimensional metric jij on each t = constant time 
slice is given by 

( 3 W = A 2 (dr 2 + r 2 dn 2 ). (2.15) 

Eqs. (2.10) and (2.11) imply 



n a = (-a, 0,0,0). 



(2.16) 



The lapse function a measures the ratio of proper time to coordinate time for a normal observer: 

dr = adt. (2.17) 

The shift /3 is the velocity of the spatial coordinates with respect to normal observers. It is a vector quantity, 
but in spherical symmetry only the radial component is nonzero (/3 a = for a ^ r), so we write /3 = /3 r . 
The shift is crucial for our AHBC method: in order for a coordinate grid point to have no causal effect on 
the region exterior to that point, the coordinate speed of radially outgoing photons must not be positive, 
and this requires a (positive) shift. 

The maximal slicing condition (2.14) is important for our conventional SA method because it causes the 
lapse function to become small in the strong-field region of a spatial slice that is about to hit a singularity. 
This slows down the passage of proper time in this region while other regions on the slice propagate farther 
into the future. In a typical black hole spacetime, maximal slices will never encounter the singularity at r = 0. 
For the AHBC method, maximal slicing is not at all necessary, but it is convenient because it eliminates a 
component of K a i, from the equations. 

We choose the isotropic spatial gauge primarily for convenience. The most general spherically symmetric 
metric can be written in the form 



ds 2 



-{a 2 - A 2 /3 2 ) dt 2 + 2A 2 /3drdt + A 2 dr 2 + B 2 r 2 dti 2 



(2.18) 



The isotropic gauge condition A = B eliminates B from the field equations. This choice is not necessary for 
either the SA or AHBC methods. However, there are gauge choices (for instance, /3 = 0) that would spoil 
the AHBC method. 

Note that Seidel and Suen[9] do not use isotropic coordinates, but instead choose the metric (2.18) and 
the gauge condition dA/dt = 0. This choice is useful because the proper distance between two coordinate 
radii remains fixed in time, but it does not seem to be necessary for the success of the AHBC method, at 
least in spherical symmetry. In fact, it tends to complicate the analysis after matching onto the SA method, 
especially in the linearized regime where it produces a frozen coordinate wave that adds gauge terms to the 
metric variables. 

In maximal isotropic gauge, the field equations (2.5) break up into two evolution equations, 

,A. 



At 

A 



P- 



a An 



(2.19) 



K T t — (3K T r = a 

and four spatial constraint equations 

1 



S_ cr / _, ^ 2 A r I A r 2 

iirb r + -A T — — + - 

4 A 3 \ A r 



^(W\ 



A 5 ' 2 



IttSL - 3K. 



8irp' + "-Kl . , 



3 

— j 

4 

(Ar), r 
Ar ' 



1 
aA 3 r 2 



(Ar 2 a : , 

'J 
r I — 

r 



-Kl + 8717/ + 4ttT', 



3 r , 
2 



2a r ( Ar ). 
A 2 Ar 



(2.20) 

(2.21) 

(2.22) 
(2.23) 
(2.24) 



Here and elsewhere in this paper, commas denote partial derivatives and 

K T = 2K% = 2K\. (2.25) 

Eqs. (2.21) and (2.22) are the Hamiltonian and momentum constraints. Eqs. (2.23) and (2.24) result from 
the maximal slicing condition K T = —K r r and the isotropic coordinate condition A = B, respectively. It 
is important to note that Eqs. (2.19)-(2.24) are not all independent: one may use the Bianchi identities to 
eliminate two of the six equations. 

The effective source terms appearing in Eqs. (2.19)-(2.24) are defined by 

p' = n a n h T' ah , 

S a — "fa n T hc , 

S' ah = 7 c a7 iZ d , (2.26) 

where T' ah is given by Eq. (2.6). Similarly, matter source terms p, T, S a , and S a i, (without the prime) are 
defined as in Eqs. (2.26), with the unprimed T a j appearing on the right hand side of each equation. If one 
replaces the primed source terms in Eqs. (2.19)-(2.24) with their unprimed counterparts, one recovers the 
Einstein equations. 



C. Scalar Field 

To solve the scalar wave equation (2.4) numerically, it is useful to define the variables 

n = -£„</>, (2.27) 

$ = <f> ir . (2.28) 

In our coordinate system, Eq. (2.27) reduces to 

^ t =/3$-oTI, (2.29) 

and the scalar wave equation (2.4) becomes 



$ t = /3$ r +/3 r $- («n), r , (2.30) 

87rTa 1 

3 + 2w ~ A 3 r 2 



n, = W,r + |^- - -^ (Ar>a*l . (2.31) 



Using II rather than <j>^ as a dynamical variable eliminates explicit time derivatives of a and /3 in 
Eq. (2.31). Furthermore, like the extrinsic curvature K a b, the quantity II is defined geometrically (by 
Eq. 2.27) on each time slice, independent of how one gets from one slice to the next: in other words, its 
value on a given slice is independent of a and /3. This is not true for <j>^. This is important when making 
the transition from the SA method to the AHBC method: the lapse and shift are discontinuous across the 
time slice that defines this transition. 

The quantity $ is a useful dynamical variable for the AHBC method because it eliminates explicit second 
derivatives from the wave equation. However, in the SA method it is sufficient to use <f> as a dynamical variable 
and to explicitly compute <f> tr and <f> trr by finite differencing when necessary. 
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Evaluating the effective source terms (2.26) using Eq. (2.6) and the metric (2.13) yields 

8 ^=(2^^( na -S)« ^ 

W* = 8,p + ^ (V + |!) + -^ ($Ar 2 ) r , (2.32b) 

8.SV* = 8^ r - ^ + ^ (n 2 + g) + HA' T + -L (|| , (2.32c) 

sti-s^ = 871-5; — - — n r - §k t . (2.32a) 



D. Linearizea Equations in Vacuum 



In the weak-fieW regime, we can use linearizea theory to aescribe the gravitational fieW. By matching 
our numerical variables to the linearizea solution, we can aetermine the gravitational raaiation seen by 
an observer at infinity, ana we can set bounaary conaitions at the outer eage of our finite-aifference gria. 
These bounaary conaitions, incluaing the ones imposea on elliptic equations, are valia even while the wave 
is passing through the bounaary. Such a matching technique was introaucea in numerical relativity by 
Abrahams ana Evans [16]. 

Define the new variables 

a = A-l, (2.33) 

e = a-l, (2.34) 

^ = <f>-l. (2.35) 

If we set all matter terms to zero in Eqs. (2.19)-(2.24) ana Eqs. (2.29)-(2.32), ana if we keep only terms 
linear in a, K T , e, /3, £, $, ana II, we obtain 

t« = ^(^,r), r , (2-36) 

n = -It, (2.37) 

*=£,r, (2.38) 

a,t = -~ \k t , (2.39) 



r 2 
2e r 2a , 



2a :r $ jr $ 
r 2 r 

3K 



$ r , (2.40) 

(2.41) 



n r , (2.42) 

(r 2 e, r ) jr = (r^ r ) ]r , (2.43) 



- r J ,r 2 

If only outgoing waves are present, the solution of Eq. (2.36) is 

r£ = f(t - r), (2.45) 
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where / is an arbitrary function of (t — r). We now insert Eq. (2.45) into Eqs. (2.39)-(2.44), and impose the 
boundary conditions e = o = at r = oo. We find the solution: 

A=l + a=l + ^-i^, (2.46) 

«=l + e=l + 4 ) + / 4 l) , (2-47) 

/3=^- / ^ ) , T (2-48) 

K T = fit ~ r) + 2fit ~ r) + ^ / (2M T + 3/(£- r) + 2C(i)) di. (2.49) 

Here C(t) is an arbitrary function of time, M T is a constant, and a prime denotes a derivative with respect 
to the argument. The gauge function C(t) results from the presence of a nonzero shift. 



E. Masses 

In order to define the mass of a spherical body in Brans-Dicke theory, it is convenient to transform the 
solution (2.45)-(2.49) into a simpler gauge. Let 

h% = h%+{ a , b + <; bia , (2.50) 

where 

I/"*/., 3 



Co=-J \M T +-f(t-r) + C(t)\dt, (2.51) 

C- = 0. (2.52) 



Here h" h l is the metric perturbation in maximal isotropic gauge, and h^J is the perturbation in Lorentz- 
Thorne (LT) gauge[16]. In LT gauge, we have 

« = ^ + ^. («3) 

KJ = 0, (2.54) 

w^m-n^l), (2 . 55 ) 



f(t - r) 
4> = 1 H -. (gauge invariant) (2.56) 

Note that there is no shift in LT gauge. 

For a static situation, / is a constant, so it will appear as an additional "mass" in the metric. We will 
therefore write 

/ = 2M S = constant, time-independent case. (2.57) 



Hence, 



*S5 = ^±^. PJM) 

I'M = ». (2-59) 

ft ,-,,(?^*), ,,.«,, 

*=1+^. (2.61, 



We see from Eq. (2.58) that a test particle in a Keplerian orbit measures a total mass M equal to 
M T + M s . The "scalar mass" M s is the portion of the active gravitational mass produced by the scalar 
field. As discussed further in Paper II, the "tensor mass" M T is the mass measured by a test black hole in 
a Keplerian orbit. In general relativity, M s = and M T = M. 

Lee[17] has derived expressions and conservation laws for the tensor and scalar masses in terms of 
superpotentials. He defines the gauge-invariant quantities 

M T =^J [^(-g) ( ff °V - ffV')]^ d 2 S,-, (2.62a) 

M T -M S =^J [(-») ( ff °y - ffV')] ,,■ d 2 S,-, (2.62b) 

which in the time-independent case reduce to M T and M T — M s , respectively. Here we assume Cartesian 
coordinates, and £; is the two-dimensional area element in the asymptotic rest frame of the source. 

Evaluating Eqs. (2.62) using the metric (2.13), we obtain 



M T =±J(<pA%r*dn = -^(<f 


' A % 


? 




(2.63a) 


M T -M a =^J(A^da = -^(AX 








(2.63b) 


In the linearized regime, we can expand these expressions to first order in 


(</,- 


1) and 


(A- 


- 1). The result is 


M T = -r 2 A :r r 2 $, 








(2.64a) 
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M s = --r J $. 

2 








(2.64b) 


Further simplification using the linearized equations (2.45)-(2.49) yields 










M T = M T , 








(2.65a) 


f rf 
M s = - + —. 








(2.65b) 



Although these reduce to M T and M s in the time-independent case, the scalar mass is not well-defined 
during dynamical epochs: M s approaches rf'(t — r)/2 as r ^ oo. The scalar mass has other unusual 
properties, which are discussed in detail by Lee[17]. For example, the scalar mass is not positive definite, 
and scalar mass carried by a gravitational wave does not curve up the background spacetime. Because the 
tensor mass M T does not suffer from such difficulties, it is the quantity that most deserves to be called 
"mass". The quantity M T is positive definite, can only decrease by the emission of gravitational radiation, 
and has other energy-like properties, unlike the scalar mass or the active gravitational mass M. 

Gauge-invariant formulas for the fluxes of M T and M s can be expressed in terms of the Landau-Lifshitz 
pseudotensor and other pseudotensors involving the scalar field [17]. In the case of spherical symmetry, it 
is easier to obtain these expressions by differentiating Eqs. (2.63) with respect to time and using the field 
equations to eliminate second derivatives of the metric. The two methods must give the same answer. The 
result, to second order in the amplitudes, is 



(2.66a) 
(2.66b) 
(2.66c) 



-2M Tjt 


= - (6A r + 3$) - -$K T + n$(w - 1) 

r ' 2 


r 2 




- (n + A' T )(4A r + a r ), 


-2(M T -M s ),t 


= (3 ( ' r $ r ) + $A T + wll$ KJAA r +a r ) 
\ r ' r ) 




+ n r (l + e-^ + 4a), 


-2A4 M 


/5$ \ 7 
= j3 1 — + $ r ] - -$A T - n$ - n(4A, r + a r ) 


r 2 




-n, r (l + t-^ + 4a), 



where the variables e, a and £ are defined by Eqs. (2.33)-(2.35), and we assume that the fluxes are evaluated 
in vacuum. Notice that M T t = to first order. This is because M T is strictly constant in linearized theory. 
However, there is a nonzero first-order contribution to M s t (the II jr term). 

If we convert the surface integrals in Eqs. (2.63) to volume integrals by Gauss' theorem, and we eliminate 
second derivatives of A using the Hamiltonian constraint, we obtain 



M T = - 
2 



o 



Sttp-M 3 + -A 6 (K T ) 2 <f> 2 + (- - l) ^A A 



--U 2 A 6 - 7(j>$A 3 A }r - 7(f> 2 A 2 (A } r ) 2 | r 2 dr, (2.67a) 



M-y-M, 



2./o 



%*~pA* 3 6 2 ^ 

—^- + - A A (A T ) + — 

$ ., . _ .,, . ,o A 4 $ r 2A 4 $ 



? 6(a ' t)2 + w {n2A2 + $2) 



-^-A 3 A t r - 7A 2 (A t r) 2 



„2 



dr. (2.67b) 



These expressions will give us an additional check on our numerical code. Because the region of integration 
contains the origin, this check is only useful in the SA method. 



F. Light Rays and Apparent Horizons 
From the metric (2.13), one can determine the coordinate velocity of radial light rays 

and the ingoing and outgoing radial null vectors 

A|=(l,±^-/3). (2.69) 

A marginally trapped surface is defined by 

V N+ A = 0, (2.70) 

where A is the area of a (t, r = const) surface: 

A = Airr 2 A 2 . (2.71) 

In maximal isotropic gauge, Eq. (2.70) takes the form 

\ + ^ = 2^- ( 2 - 72 ) 

where we have used Eq. (2.19) to eliminate the time derivative of A. 

The apparent horizon is the outermost marginally trapped surface. In our standard SA method, we use 
Eq. (2.72) to locate trapped surfaces. In the AHBC method, we place a grid point on the apparent horizon 
and use Eq. (2.72) as a boundary condition at that grid point. 

For the AHBC method, we also need a relation that forces the apparent horizon to remain at a constant 
grid point as the metric variables and scalar field evolve in time. Such an equation can be obtained by 
differentiating Eq. (2.72) with respect to t (holding r constant), and then using Eqs. (2.19)-(2.22) to eliminate 
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time derivatives, second radial derivatives of A, and first radial derivatives of K T . After some algebra, we 
arrive at the result 

a_ = l-8irA 2 r 2 (p' + S' r /A) 
PA l + 8TrA 2 r 2 (S r r i + S> r /AY [ ' ' 

where the primed source terms are given by Eq. (2.32). 

III. SINGULARITY-AVOIDING (SA) METHOD 



Here we present our singularity-avoiding numerical method for solving the Brans-Dicke equations for a 
spherically symmetric distribution of dust. Except for the addition of the Brans-Dicke source terms in the 
field equations, this method is very similar to that of Ref. [18]. 

We use a mean field particle simulation scheme in which the collisionless matter distribution is sampled 
by a finite number of noninteracting particles, and the gravitational field variables are defined on a finite 
number of grid zones. At each time step, the particles are moved according to the geodesic equation. By 
binning the particles, we calculate the source terms p, S r , and T in each grid zone. The source term T is 
used to determine the scalar field via the wave equation, which is solved by finite differencing on the grid. 
All three source terms appear in the constraint equations, from which we determine the metric variables. 
The process is then repeated for the next time step. 



A. Matter Source Terms 

We sample the matter by a finite number of noninteracting particles. If each particle has rest mass m, 
4- velocity u a , and comoving number density A/", the stress-energy tensor is given by 

T ah =^2mAfu a u h , (3.1) 

where the sum is over all of the particles. The source terms defined by Eq. (2.26) are 

p= ^2mAf(au ) 2 , (3.2) 

T=-^mA', (3.3) 

S r = -^2mAf(<xu°)u r , (3.4) 

S r r = ^ j mtiu 2 lA 2 . (3.5) 

The comoving number density A/" of each particle is proportional to a delta function for point particles, 
but can be treated as a continuous quantity if the number of particles is large and the particle distribution 
is smooth. In evaluating the source terms, we average over all particles in each grid zone, so that we can 
think of each particle as occupying the entire volume of the zone. Therefore, we write 

= 47r(aM°)A 3 r 2 Ar, (3.6) 

where Ar is the width of a grid zone. 

Notice that the source terms in Eqs. (3.2)-(3.5) depend on the metric function A. Since we use these 
source terms to solve for the metric functions, it is helpful to define quantities that can be calculated from 
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the particle variables alone, without referring to the metric. This is especially important when solving the 
nonlinear Hamiltonian constraint (2.21) for A. We define 

p = pA* = Y m (°f\ (3.7) 

T = TA 3 = -J2 4ir(au0)r2Ar , (3-8) 

L^i 47r r 2Ar v ' 

~S\ = S r r A 5 = Y , mU } „ . (3.10) 

Although the metric function a still appears in these expressions, it is always multiplied by u°. The quantity 
au° depends only weakly on the metric (Eq. 3.12). 

Simply evaluating the sums in Eqs. (3.7)-(3.10) in each grid zone gives reasonable values for the source 
terms in that zone. Even better is to use the algorithm of Hockney and Eastwood[19] to distribute each 
particle's rest mass among its own grid zone and each of the two neighboring zones. This procedure reduces 
the stochastic fluctuations associated with having only a finite number of particles, giving us smoother 
results[20]. 



B. Geodesic Equations 



Since the matter is made up of particles, the equation of motion V a T ai = reduces to the geodesic 
equation for each particle. Using the metric (2.13), we write the geodesic equation as 

dr au r 

— = — — — — — a (3.11a) 

-=>M«, r + ^ + ^^+ (M0 ^(- + T j. (3-llb) 

The normalization of four- velocity gives us 

These equations are solved for the variables r, u r , and au° for each particle at each time step by an 
embedded fourth-fifth order Runge-Kutta scheme with adaptive step size[21]. The quantity u^ is a constant 
of the motion. Derivatives of the metric that appear in Eqs. (3.11)— (3.12) are obtained on the grid by 
finite differencing. The values of these derivatives and the metric functions at the particle position are then 
determined by interpolating from the grid. 



C. Scalar Wave Equation 

After moving the particles and determining the source terms, we proceed to the wave equation. In order 
to minimize roundoff error in the nearly GR regime ((f) close to unity), we use the variable 

£ = ^-1. (3.13) 
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The wave equation can then be written 

£,* - 2/3r^ r 2 = -an, (3.14a) 

n t -2/3rn r2 = 8lTTa --fL-(r 3 Aa£ r3 ) , . (3.14b) 

^ ' r A 3 (3 + 2w) A 3 V Sr V 3 K ' 

These are Eqs. (2.29) and (2.31) written without explicit use of the variable $. It is not necessary to use $ 
as a dynamical variable until we introduce the AHBC method in Section IV. 

The regularity condition at the origin is 

£ ir =n r = 0. (3.15) 

Since £ and II behave like 

C 1 + C 2 r 2 , C 1 ,C 2 constant (3.16) 

for small r, we take radial derivatives with respect to r 2 rather than r in Eqs. (3.14). This ensures that our 
finite difference equations yield the correct result near the origin [22]. 

The outgoing-wave condition at infinity is 

(rY\ r + (rY\ t = 0, (3.17) 

where Y is either II or £. For spherical symmetry, Eq. (3.17) is exact in linearized theory, as one can verify 
by substituting Eq. (2.45). 

We solve Eqs. (3.14) by an explicit staggered leapfrog scheme that determines £ and II at timestep n + 1 
given their values at timesteps n and n — l. No information at timestep n + 1 is needed to solve the equations. 
After determining £ and II, we calculate <f> from Eq. (3.13). Finite difference approximations are presented 
in Appendix A.l. 

D. Constraints 



After determining the quantities £, II, and <f>, we then solve for A and K T . 
By introducing the new variables 

Z = A 3 r 3 <l)K T , (3.18) 

V-^A 1 / 2 , (3.19) 

we can rewrite the momentum constraint (2.22) and the Hamiltonian constraint (2.21) in the form 

871-5; , fi 2t/> 6 wII£,.2 

5Z r 5 = - 2V> 6 II r2 - — — -^, (3.20) 

' r ' 6 



6 'rV'YrO =-tAw1S# 2 ^ 



-^f -^f%^V--(^rO.- (3-21) 

In general relativity II = £ jr = 0, so one can determine Z from Eq. (3.20) and then compute t/> from Eq. (3.21). 
This is the reason for using the variable Z instead of K T . In Brans-Dicke theory, Eqs. (3.20) and (3.21) are 
coupled because of the last two terms in Eq. (3.20), so we must solve these equations simultaneously. 

13 



Eq. (3.20) requires a single boundary condition. For regularity, we set Z = at the origin. Eq. (3.21) 
requires two boundary conditions. The regularity condition at the origin is 

V> = 0. (3.22) 

At infinity, we match to the linearized solution (2.46). Differentiating this equation yields 

(r^), r = l-^, (3.23) 

which does not require knowledge of the tensor mass M T . 

We determine the variable Z using the momentum constraint (3.20), and we solve the nonlinear Hamil- 
tonian constraint (3.21) for t/> using the iterative scheme described in Appendix A. 2. Because t/> appears in 
Eq. (3.20), we must recompute Z at each step in the iteration. We obtain an initial guess for t/> from the 
evolution equation (2.19), which we write as 

V>,t - /3V> = ^ - ^«VAV (3-24) 

After determining Z and t/>, the variables A and K T are found from Eqs. (3.18) and (3.19). Finite difference 
forms of Eqs. (3.20)-(3.24) are presented in Appendix A. 2. 

E. Lapse and Shift 



Having determined the scalar field variables and the spatial metric variables, we can now calculate a 
and /3. The lapse equation (2.23) can be written 

6/i, x 6 , „ „ N 8tt / f \ wll 2 3 T ^ . 

^ ( rAa ^i> = W , O^O,,- + J [<> + 2T37^ J + !jr + 2 A - (3 - 25) 

The regularity condition at the origin is 

a r = 0. (3.26) 

We match to the linearized solution to obtain the boundary condition at infinity: by differentiating Eq. (2.47), 
we obtain 

(ra), r = 1 + rlL (3.27) 

This condition is independent of the gauge function C(i) that appears in Eq. (2.47). 

After solving for a, we calculate /3 from the shift equation (2.24). This equation requires a single 
boundary condition. We impose Eq. (2.48) at the outermost grid point: 

' 2 2 K J 

Finite difference forms of Eqs. (3.25)- (3.28) are given in Appendix A. 3. 



F. Grid 



Our numerical grid extends from r\ = to r 8max = r max . All variables are centered at half-grid points 
r i+i/2- In order to obtain a nearly constant number of particles in each grid zone, we divide the grid into 
inner and outer regions. The inner region, which contains all the particles, extends from r\ to r 8 - art where 
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r-{ art is the grid point just outside the outermost particle. The grid point r-2 is chosen so that it contains a 
fraction l/i pa rt of the rest mass. Spacing between other grid points in the inner region is geometric in r 3 : 



i + l '% _ ' i + 2 






(3.29) 



This ensures that the grid spacing is close to uniform in r 3 , so that for a uniform particle distribution each 
zone has approximately the same number of particles. The outer grid, which extends from r 8 - art to r 8max , is 
matched smoothly onto the inner grid. The spacing between outer grid points is also geometric in r 3 . 

The grid is allowed to move at every time step so that it follows the particle distribution. Each time the 
grid is moved, all gravitational field variables are interpolated from the old grid onto the new one. Moving 
the grid allows us to place many grid points where they are required to maintain accuracy. If the grid 
remained stationary during gravitational collapse, the particles would soon end up only in the innermost 
grid zones, invalidating our finite difference approximations. 



G. Identifying Apparent Horizons 



To determine the location of an apparent horizon, we evaluate the quantity 

9 = i + 2%-iv 2 ^ T (3.30) 

r ip 2 

at each grid point. By Eq. (2.72), all grid points with < are contained within a trapped surface. 
Therefore, if r 8AH is the outermost grid point with 6 < 0, then the apparent horizon lies between r 8AH and 
r 8AH _|_i. We obtain the approximate value of tah by 

r AH = r iAH - 6 iAH - — . (3.31) 

" AH + l "'AH 



H. Initial Data 



Our initial time slice occurs at a moment of time symmetry, so that K T = II = /3 = 0. We first calculate 
the initial values of the metric and scalar field using a method similar to that of Matsuda[23]. This method 
is discussed in Appendix B. It involves only ODEs, which can be solved to arbitrary accuracy by standard 
numerical methods. 

After obtaining the ODE solution, we re-solve for the initial data using the mean-field particle scheme. 
We first randomly place particles in the interior according to the rest mass function M Test (r) obtained from 
the ODE solution. We then compute the matter source terms by binning the particles, and we solve for ^ 
using Eq. (3.14b ) with II = /3 = 0. Next we solve for t/> and a using the Hamiltonian constraint (3.21) and 
the lapse equation (3.25). By comparing £, t/>, and a with the result of the ODE solution, we obtain an 
important check on our method. 



I. Diagnostics 

There are two evolution equations, Eqs. (2.19) and (2.20), that we did not use in solving the Brans-Dicke 
equations (Although we used Eq. (2.19) as an initial guess for t/>, we refined our guess using the Hamiltonian 
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constraint). We can test the accuracy of our code by comparing the left and right-hand sides of these 
equations at each time step. 

Another check is obtained by mass conservation. In the linearized region far from the origin, we calculate 
the quantities M s and M T by Eqs. (2.65), and compare with the result given by Eqs. (2.66) and with the 
volume integrals (2.67). This is more useful than evaluating the evolution equations at each time step because 
it is sensitive to errors that accumulate over time. 



IV. APPARENT HORIZON BOUNDARY CONDITION (AHBC) METHOD 



Here we present our AHBC method of solving the Brans-Dicke equations for a spacetime that contains 
a black hole. Our method depends on two primary ingredients. The first is the existence of what we call 
a coordinate causal horizon (CCH), which is a coordinate surface inside of which the coordinate velocity of 
radially outgoing light rays is negative, so that instantaneously any point inside the CCH cannot causally 
influence a point in the exterior. Note that a CCH is coordinate dependent, and can exist even in Minkowski 
space by choosing a coordinate system that moves faster than light. 

The second ingredient is a horizon-locked coordinate system, similar to that of Seidel and Suen[9]. We 
place the AH at a fixed radial coordinate r^n, and require it to remain there at all times. 

We truncate our spacetime at some radial coordinate r± such that < r± < r^n, and we discard the 
interior region. As long as a CCH exists for some r > r±, discarding the interior cannot in principle affect 
the evolution of the exterior. In general relativity, a horizon-locked coordinate system guarantees that the 
AH is on or inside a CCH, since a light ray at the AH cannot propagate outwards. In this case, one can 
set r*i = tah- In Brans-Dicke theory, we find that this is not the case, as discussed further in Paper II. 
We therefore leave a small buffer zone between r± and r^n, and check that a CCH is always present in the 
portion of spacetime that we retain. Even in general relativity, this buffer zone is also useful for computing 
radial derivatives at the AH. 

We solve the wave equation in a manner that takes advantage of the causal structure of the spacetime. 
We define a "causal boundary" at some radial grid point tqb which coincides with either the AH or the 
CCH, whichever is smaller. Given any r' < tqb, our difference scheme ensures that the scalar field at r' 
depends only on quantities at r > r' . Furthermore, no explicit boundary condition is needed at the inner 
grid point r±. 

To obtain the metric variables a, A, /3, and K T , we solve spatial constraint equations that require a 
total of six boundary conditions. Three of these are provided by matching to the linearized solution at 
the outer grid point. The other three, in the SA method described in the last section, were obtained by 
regularity at the origin. In the AHBC method, however, we exclude the origin from our spacetime, so these 
three boundary conditions must be imposed at the AH. We obtain two conditions by locking the AH to a 
fixed radial grid point. One is the marginally trapped surface relation (2.72), and the other, Eq. (2.73), is 
the requirement that the AH remain at a constant coordinate radius for all time. A third condition could 
be obtained by setting the tensor mass, M T , of the black hole to the value obtained by mass conservation, 
Eq. (2.66a ). However, this would prevent us from using mass conservation as a check on our numerical 
integrations. Instead, we use a different approach: we use the evolution equation (2.19) to update A on the 
AH. Another possibility would be to evolve the value of K T on the AH using Eq. (2.20), but this would be 
more complicated. 



A. Matter 



As in the SA method, the particles can be moved according to the geodesic equations and the matter 
source terms can be computed by binning particles into zones. However, in the case of gravitational collapse 
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starting with a uniform density profile, we will find that an AH does not form until after all matter has 
fallen into the black hole. For this reason, we only need to solve the AHBC equations in vacuum — it is not 
necessary to include matter particles in the code. This is a great advantage in terms of efficiency, since in 
the SA scheme most of the computer time is spent moving the particles, and it also gives us flexibility in the 
grid choice, since we no longer need to choose grid zones that are approximately constant in volume. 

We emphasize that there is no fundamental restriction that prevents us from including matter in the 
AHBC method, and for some physical situations not investigated in this paper, e. g., the collapse of a 
distribution with a core-halo structure, an AH may form while there is still matter in the exterior region. 
In such a case, including particles in the code would not be very difficult. For this reason, in the following 
sections we will continue to include the matter source terms in our equations. 



B. Grid 



For accuracy, it is usually desirable to choose the finite difference grid to be uniformly spaced. However, 
a grid uniformly spaced in r does not yield enough grid coverage near the AH. This is because the initial data 
for the AHBC method is determined by matching onto the SA method after a horizon is formed, and the SA 
method typically produces an apparent horizon radius tah much less than the outer grid radius and metric 
quantities that vary exponentially with r near the AH. It is therefore convenient to choose a logarithmically 
spaced grid, and to take all radial derivatives with respect to In r rather than r. This choice places many 
grid points near the AH where they are needed, and still allows us to take finite differences on a uniformly 
spaced (in In r) grid. 

For convenience, we define the new variable 

r\ = In r. (4-1) 

Our radial grid then consists of the i max points (r/i, r/2, . . . , f]i AH , . . . , ?ft max ). The AH is located at the 
grid point i = i'ah- The point r 8max is chosen to be far from the black hole, in a region where linearized 
theory is valid. All variables except $ are defined on the grid points; $ is defined on the half-grid points 
(r?3/ 2 , • • • , % max +i/2)- The staggering of $ is essential for our method of solving the wave equation. 

Unfortunately, spacing the grid uniformly in In r tends to produce too sparse a grid at large radii. A 
consequence of this is that outgoing waves get partially reflected at the outer edge of the grid. To minimize 
this effect, we pick a maximum threshold Ar*o for the spacing between grid points. We choose our grid 
uniform in r\ unless this choice would yield a grid spacing (in r) greater than Ar*o at the outer boundary. In 
the latter case, we set r 8max — r 8max _i equal to Ar*o, and choose the other grid points so that 

Vi+i- Vi _ Vi+2 ~ Vi+i / 4 2 x 

r?i-r?i_i r] i+1 -r]i 

This gives us a grid that is as close as possible to being uniform in r\ while still obeying r 8 _|_i — r 8 - < Ar*o for 
all i. 



C. Wave equation 



The wave equation can be written in the form 

= /3$-«n, (4.3) 

*t = -[/3$, +/3,,$-(«ny, (4.4) 



r 



n, t = -n 



/3 87rTa a 
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,v 



a ip 



(4.5) 



We solve Eqs. (4.4) and (4.5) for II and $ by the causal method described in Appendix C. Because this 
method takes into account the causal structure of the spacetime, it is not necessary to impose explicit 
boundary conditions at the inner edge of the grid. Boundary conditions at the outer edge of the grid are 
obtained by matching both II and $ to the outgoing-wave linearized solution. Eqs. (2.37), (2.38), and (2.45) 
yield 

1 2$ 

n t + -$,„ + — = 0, (4.6) 

i 2$ n 

$ >* + ^ + v"7 = ' (4 ' 7) 

After determining $ and II, we solve for £ using Eq. (4.3). This equation is an ODE in time and requires 
no boundary conditions. The scalar field <f> is then calculated from Eq. (3.13). 



D. Constraints 

After solving the wave equation, the Hamiltonian and momentum constraints are solved simultaneously 
for the variables t/> and Z, from which we can compute A and K T . These equations are coupled because 
scalar wave terms containing t/> appear in the momentum constraint (2.22). Eqs. (2.21) and (2.22) can be 
written as 

Z tV =87rr 4 5 r - V 6 r 3 (n,„ + 5^ , (4.8) 

$r\ 3 Z 2 2v~pr 2 ujip 5 Tl 2 r 2 



i>, m + V 1 ,!) 1 



2(f) J I6(f> 2 ip 7 r 4 (f)ip 

uiip<f> 2 r 2 rip 

8^2 4^ 



($,,+2$). (4.9) 



Eq. (4.8) requires a single boundary condition. At the apparent horizon, we impose the marginally 
trapped surface condition (2.72), which we rewrite using the variables t/>, Z, and r/: 

^ + I = 4VW (4 - 10) 

Eq. (4.9) requires two boundary conditions. At the outermost grid point, we match to linearized theory 
by imposing Eq. (3.23): 

r 2 n 

(ri>),r, = r- — . (4.11) 

Because there is no extra equation for the inner boundary condition, we use the evolution equation (2.19) 
to compute the value of t/> at the apparent horizon. This equation takes the form 

^ = A + A _ I a * T . (4.12) 

t/> ipr 2r 4 v ' 

Eqs. (4.8) — (4.11) comprise a coupled system of nonlinear equations. We solve them simultaneously by 
the iteration scheme described in Appendix D. 



E. Lapse and Shift 



After obtaining t/> and Z, we solve the lapse and shift equations simultaneously for a and /3. In the SA 
method, these equations are solved independently; here they are coupled because of a boundary condition 
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that we now impose at the apparent horizon. We write these equations as 



,7777 \ iV 






P 



3 (v ,2A2^87rr 2 I p + T 
- (K T rip z ) + 



uip 4 Jl 2 r 2 



-aKn 




At the outer grid point, we impose Eqs. (3.27) and (2.48), which take the form 



(ra\, 



P 



r + r n, 



2 2 

For the inner boundary condition on a we use Eq. (2.73), which can be written as 

a _ 1/r 2 - F x - F 3 



/3t/> 2 1/^ + ^2 + ^3 + ^4' 



where 



Fi 



Fo 




cr 



F 4 






f 3 = ^ (nv> 2 - $ 



3 + 2w V 2 / ' 

2 ^ 2 n,„ | $ 

r<j) r<j) r<j>i]) 



v _l -*,n _ 2 ®i ) ,v 



(All - $) . 



(4.13) 
(4.14) 

(4.15) 
(4.16) 

(4.17) 

(4.18a) 

(4.18b) 

(4.18c) 
(4.18d) 



This condition forces a grid point to remain at the apparent horizon. It also couples the lapse and shift 
equations (4.13) and (4.14), so that we must solve these equations simultaneously. Finite-difference forms of 
Eqs. (4.14)-(4.18) are presented in Appendix D. 

Notice that in vacuum with constant scalar field, Eq. (4.17) reduces to a = A/3. Combining this with 
Eqs. (4.10) and (4.12), we see that in this case t/> jt = on the apparent horizon. As a result, all quantities 
are manifestly time independent: Given the value of t/> on the AH, which is constant in time, the quantities 
A, a, /3, and K T are completely determined by coupled spatial constraint equations with no source terms. 

This time independence should come as no great surprise because the vacuum solution with constant 
scalar field is simply the exterior Schwarzschild solution. However, this result is remarkable from a practical 
point of view: When one uses standard singularity-avoiding schemes (including our SA scheme) to compute 
spherical collapse in general relativity, one does not obtain a time-independent system at the end of the 
simulation. Instead, one finds that the metric functions and their derivatives are changing exponentially 
with time inside the AH. However, using the AHBC method, any system that results in a Schwarzschild 
black hole will become manifestly time independent. This is the great benefit of the AHBC method: One can 
integrate black-hole spacetimes arbitrarily far into the future without encountering singularities and without 
causing metric functions or their derivatives to blow up. 



F. Initial Data 



The initial time slice for the AHBC method can be any maximal isotropic time slice that contains an 
apparent horizon. We use a slice provided by the SA method. Because <f> is a scalar and t/>, K T , and n are 
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three-dimensional geometric objects on the slice, and because we use the same 3-metric for both the SA and 
AHBC methods, the values of these variables on the initial slice can be simply read off from the SA slice. 
After laying down a new spatial grid as described in section A, we spatially interpolate these variables from 
the old SA grid to the new one. We then determine $ from Eq. (2.28). For self-consistency, we then freeze 
the value of t/> at the AH and recalculate t/> and K T from Eqs. (4.8)-(4.11). Finally, we calculate a and /3 by 
solving Eqs. (4.14)-(4.18). 

Although the initial AHBC time slice coincides with an SA slice, subsequent AHBC slices do not coincide 
with those that would have been produced by continuing the SA method, because the two methods have a 
different lapse and shift. These differences are due soley to the different inner boundary conditions used in 
the two methods. 



G. Diagnostics 

In determining the metric and the scalar field, we never use the K T evolution equation (2.20) or the 
definition of $, Eq. (2.28). Furthermore, while we need the A evolution equation (2.19) to determine a 
boundary condition at the AH, this equation is not used at any other grid point except as an initial guess. 
We therefore have three equations that can serve as diagnostics at each time step. 

In addition, we can compare the results of the masses calculated by Eqs. (2.65) and (2.66). Unlike 
the step-by-step comparison of equations (2.20), (2.28), and (2.19), this method is sensitive to errors that 
accumulate over time. Because we exclude the origin from our spacetime, we cannot evaluate the volume 
integrals (2.67) in the AHBC method. 

V. NUMERICAL RESULTS 



In this section we treat a few select collapse scenarios to calibrate our code. We use these cases to 
compare the standard SA scheme with the AHBC method. For a more complete summary of physical results 
and an evaluation of collapse to black holes in Brans-Dicke theory, see Paper II. 

A. Oppenheimer-Snyder Collapse in General Relativity 

A useful test of our code is Oppenheimer-Snyder[24] collapse to a black hole in general relativity. We 
start with a momentarily static, uniform sphere of dust with an initial areal radius R s = 10M, where M 
is the total mass of the configuration as measured by Keplerian orbits far from the origin. We follow the 
collapse of this object to a black hole using our SA scheme with w = 10 37 . Such a large value of w puts 
us well into the general relativistic limit. We use 41 interior grid points, 87 exterior grid points, and 1200 
particles. Our outer grid boundary is at r = 100M. 

Figures 1, 2, and 3 show snapshots of metric parameters at selected values of coordinate time. Also 
shown are the exact results obtained[25] by transforming the analytic solution[24] into the maximal isotropic 
gauge. The SA method agrees well with the exact solution until very late times. The "collapse of the lapse" 
can be seen in Figure 1: inside r s = 1.5M, the lapse function decreases exponentially with time, freezing 
the clocks of normal observers in this region. The values of A and /3 inside r s = 1.5M also change rapidly 
with time, as can be seen in Figures 2 and 3. Because of the large values of A, coordinate radii become very 
small compared to isotropic radii. For example, the areal radius of the matter surface at t = 87M is about 
r s = 1.5M, while its coordinate radius is only r ~ 5 x 10 _8 M. In addition, the gradients of a, A, and /3 
become very steep near r s = 1.5M. 

Tracking particles at r ~ 10 _8 M with a grid that extends to r = 100M and resolving the metric function 
A which drops from > 10 7 to 2 between r = 10 _7 M and r = M requires an enormous dynamic range. 
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Figure 1 The lapse function a versus areal radius r s on selected maximal time 
slices for general relativistic Oppenheimer-Snyder collapse from R s = 10M. Cal- 
culations by the SA method (dotted line) are compared with the exact solution 
(dashed line). Time and distance are measured in units of M . 
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Figure 2 The metric component A 1 ' 2 versus areal radius r s on selected maximal 
time slices. Labeling is the same as in Figure 1. 



Because our inner grid follows the particles at each time step, grid points tend to accumulate at the limit 
(matter) surface, r s ~ 1.5M. Although this allows us to resolve the large metric gradients better than with 
a stationary grid, our code eventually loses accuracy as the gradients grow. This can be seen in Figures 1-3: 
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Figure 3 The radial shift /3 versus areal radius r s on selected maximal time slices. 
Labeling is the same as in Figure 1. 




Figure 4 Oppenheimer-Snyder collapse in maximal slicing, calculated by the SA 
method. The five solid lines represent world lines of matter elements containing 
20%, 40%, 60%, 80%, and 100% of the interior rest mass. Time and areal radius 
are measured in units of M. The dotted line is the apparent horizon, which forms 
at about t = 44M, at which point it already coincides with the event horizon. 



the SA method and the exact soultion begin to disagree at about t = 60M, and this discrepancy increases 
exponentially with time. Even if we were able to maintain accuracy, numerical overflow would eventually 
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cause our SA code to terminate. 

A spacetime diagram of Oppenheimer-Snyder collapse is shown in Figure 4. Because of the "collapse of 
the lapse", the matter particles sit at a constant areal radius after t = 50M. An apparent horizon, which 
forms at t = 44M outside of the matter, has an initial areal radius of r s = 2M , in agreement with the 
Schwarzschild solution. It increases slightly in size after about t = 70M because of numerical errors. 
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Figure 5 Metric quantities for t > 45M for Oppenheimer-Snyder collapse in gen- 
eral relativity, calculated by the AHBC method. We remove the region of spacetime 
inside the apparent horizon at r s = 2M and only retain the vacuum exterior. All 
quantities are constant in time. 



The above numerical difficulties do not occur in the AHBC method because we no longer use such a 
pathological coordinate system. In Figure 5 we show the metric coefficients obtained from the AHBC method 
after matching onto an SA time slice at t = 45. We use 128 grid points, and place the outer boundary at 
r = 100M. The apparent horizon is located at the grid point i = 4, which is at coordinate radius r = 0.78M 
and areal radius r s = 2M. The innermost grid point is at r = 0.7Mo. Because all the matter has fallen past 
the AH by t = 45M, we no longer need to move the particles — we have discarded the region of spacetime 
that contains them. In this case the metric is manifestly static. 



B. Scalar Waves on a Schwarzschild Background 



Consider a small scalar perturbation about a Schwarzschild black hole. If the perturbation is so small 
that its contribution to the metric is negligible, then it is not necessary to solve the coupled Brans-Dicke 
equations (2.4) and (2.5) to determine the future evolution of <f>. Instead, one only needs to solve the wave 
equation (2.4) in vacuum on a Schwarzschild background. This equation can be written in the form[26, 27] 



Vu, 



(5.1) 



where 
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u = r s <f>, (5.2) 

z = r s + 21n(r s /2-l), (5.4) 

r s is the areal radius, and t s is the Schwarzschild time coordinate. The variable z is the familiar "tortoise" 
coordinate, which runs from z = — oo (at r s = 2) to z = oo (at r s = oo). The mass of the black hole has 
been set equal to unity. Eq. (5.1) provides an important check for our AHBC code. This equation is not 
difficult to solve numerically — it is simply a one-dimensional flat space wave equation in a static potential. 

We set up the following test case: At t s = 0, the metric is Schwarzschild and the scalar field is given by 

(5.5) 

where r*o = 80, C = 10 -6 , and a = 5. We set <j)^ s = initially. As time progresses, the initial pulse divides 
into two pieces: one moves outwards to infinity, while the other moves toward the black hole and is partially 
reflected by the Schwarzschild potential. 

We calculate <f> as a function of r s and t s for this case by two independent methods: we solve the 
full Brans-Dicke equations using the AHBC method, and we solve the perturbation equation (5.1) using a 
staggered leapfrog finite difference scheme. The AHBC method requires metric coefficients on an initial time 
slice; these are provided by matching onto Oppenheimer-Snyder collapse (calculated from the exact solution 
following Petrich et al. [25]) after the matter has fallen into the black hole. 

Because the AHBC method and the perturbation method use different coordinate systems, we must be 
careful in setting up the initial data. Eq. (5.5) is defined on the Schwarzschild time slice t s = 0. This slice 
does not correspond to the maximal slice t = 0, or any other t = constant slice in maximal isotropic gauge. 
However, the initial Gaussian pulse (5.5) is far from the black hole, where t = constant and t s = constant 
slices nearly coincide. Therefore, if we choose the slice t = to coincide with t s = at r s = 80, we can use 
Eq. (5.5) to determine <f> at t = for the AHBC method without introducing much error. 

Although t = constant and t s = constant slices nearly coincide at r s ~ 80 where the initial wave 
is nonzero, these slices differ considerably in the strong field region near the black hole. If we wish to 
compare the results of the AHBC and perturbation methods in this region, we must use coordinate invariant 
quantities. We therefore introduce a set of stationary observers at specified values of the areal radius r s , and 
we record the values of <f> and <f> tT measured by each of these observers as a function of r, the proper time 
measured by the clock of each observer. To synchronize the clocks in a manner independent of the choice of 
time slicing, we introduce an ingoing light signal that passes r s = 80 at t s = t = 0. Each observer starts his 
clock running from r = when he sees the signal, e. g., at the initial time slice, an observer at r s = 80 reads 
r = 0, an observer at r s = 90 reads r ~ 10 (plus 0(l/r s ) corrections), and an observer at r s = 70 reads 

T 10. 

Figure 6 shows <f> and <f> tT as measured by two different observers. The observer at r s = 100M sees two 
peaks of scalar radiation. The first travels outward from r s = 80M, starting at the observer's proper time 
r ~ 20M, and passes the observer at r ~ 40M. The second travels inward from r s = 80M, is partially 
reflected by the black hole, and then moves outward, passing the observer at r ~ 220M. The observer at 
r s = 5M does not see the outgoing radiation pulse, but instead sees a combination of the ingoing pulse 
and its reflection. The agreement between the two numerical methods is excellent, demonstrating that the 
AHBC scheme is able to handle gravitational radiation without producing large numerical reflections at the 
apparent horizon or at the outer grid boundary. In addition, the above case extends from t = to t = 300M, 
much longer than a traditional SA scheme would allow. 

Figure 7 shows the coordinate velocity of outgoing light rays, 

5 = 1-* <-» 
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Figure 6 Scalar field and its derivative with respect to proper time r, as measured 
by stationary observers, for a small scalar perturbation about a Schwarzschild black 
hole. The proper time r is measured in units of M, the mass of the black hole. 
The AHBC method (solid line) and the perturbation method (dotted line) almost 
coincide. Both methods use 256 grid points. The observation radii, from top to 
bottom, are r s = 100M and 5M. 



This quantity is negative inside the apparent horizon, indicating that in this region information can only 
move inward with respect to the coordinates. Such a coordinate system is essential for solving the wave 
equation without explicitly imposing a boundary condition at the inner edge of the numerical grid. Because 
the spacetime is Schwarzschild, the apparent horizon coincides with the event horizon, so that the coordinate 
grid point on the AH moves along an outgoing light ray. At large radii dr/dt approaches unity because the 
coordinate system is asymptotically Minkowskian. 



C. Oscillating Einstein Cluster 



In order to demonstrate the emission of monopole gravitational radiation, we evolve a system that 
undergoes spherical oscillations. At t = we place particles in randomly oriented stable circular orbits, and 
then reduce each particle's four-velocity component ug by a constant factor £. The particles are arranged 
so that the comoving energy density is initially uniform throughout the configuration. If £ = 1, the particle 
distribution would remain in equilibrium. In general relativity such an equilibrium system is called an 
Einstein[28] cluster. For £ < 1 and a weak gravitational field, particles do not remain in equilibrium, but 
instead move on elliptical orbits with identical periods, and the entire spherical distribution periodically 
expands and contracts in a homologous manner[29]. In strong gravitational fields, however, the oscillations 
of the distribution are not homologous because of shell crossing. In addition, the existence of a scalar field 
in Brans-Dicke theory allows monopole radiation, so that the system loses energy. The particles eventually 
settle into a new equilibrium state. 

Figure 8 shows <f> — 1, II, and the metric quantities A 1 ' 2 — 1 and a — 1, measured at r = 400Mo as a 
function of coordinate time for an oscillating cluster with £ = 0.95 and an initial areal radius R = lOMo. 
Here Mq is the initial active gravitational mass M of the cluster. The particles and the spacetime are evolved 
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Figure 7 Coordinate velocity dr/dt of outgoing light rays versus areal radius 
r s for a Schwarzschild black hole with a weak scalar perturbation. Areal radius 
is measured in units of M, the mass of the black hole. The quantity dr/dt is 
time-independent because the scalar field is too small to significantly change the 
background Schwarzschild metric. The arrow indicates the position of the apparent 
horizon, which is at r s = 2M. The inner edge of the computational domain is at 
r, = 1.96M. 



using the SA method. In order to produce a significant amount of gravitational radiation, we set to = 1. 
As the cluster oscillates, it produces scalar gravitational radiation that first reaches r = 400Mo at a time 
t ~ 400Mo. The amplitude of the oscillation is damped because energy is lost to radiation. Notice that the 
lapse function a oscillates even before the radiation reaches r = 400Mo. This is because a is determined by 
elliptic equations which "feel" the oscillations of the interior matter distribution as well as the gravitational 
radiation that propagates to infinity. 



D. Oppenheimer-Snyder Collapse in Brans-Dicke Theory 



We now consider a case in which the effects of the scalar field are large enough that one cannot use linear 
approximations to determine the evolution of the spacetime. We treat the analogue of Oppenheimer-Snyder 
collapse in Brans-Dicke theory. This case is identical to the one in section A except that we set w = 1 to 
increase the effect of the scalar field on the spacetime geometry. This value of to conflicts with experiment; 
however, our purpose here is not to calculate a physically realistic collapse, but to provide a test of our 
numerical method. 

We begin with a uniform, momentarily static sphere of noninteracting particles with total mass Mq and 
initial radius r s = lOMo. The zero subscript on the mass refers to the initial value — as it evolves in time, the 
system loses mass because of scalar gravitational radiation. Because to = 1, the scalar mass M s contributes 
about 16% of the initial mass Mo. We follow the collapse of the configuration with our SA method, using 81 
interior grid points, 175 exterior grid points, and 1200 particles. Our outer grid boundary is at r = lOOMo. 

A spacetime diagram of the collapse is shown in Figure 9. The apparent horizon forms outside the 
matter surface at t = 44Mq. It initially has a radius r s = 1.42Mq, but it expands as the scalar field radiates 
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Figure 8 Scalar field and metric quantities at r = 400Mo for an oscillating 
Einstein cluster with to = 1, £ = 0.95, and an initial areal radius R = lOMo. 
The solid lines show results calculated with an outer grid radius at r = 410Mo, 
1200 particles, 61 inner grid points, and 195 outer grid points. The dotted lines 
show results calculated with an outer grid radius at r = lOOOMo, 1200 particles, 61 
inner grid points, and 451 outer grid points. The close agreement between the two 
cases demonstrates that our code accurately handles the outgoing wave boundary 
conditions at our outermost grid point, even while waves are passing through the 
boundary. 



away. By t = 65Mo the AH has grown to r s = 1.69Mo It then decreases by about O.OlMo (not visible in 
the figure) and reaches equilibrium. This decrease represents a feature of Brans-Dicke black holes that is 
not present in general relativity. It is discussed further in Paper II. After t = 90Mo, numerical errors from 
exponentially growing metric gradients become large enough that the matter and AH both move outward 
in areal radius. These errors eventually cause the code to crash at t = 99Mo. Adding more grid points does 
not allow the code to run much longer: with twice as many grid points, it only runs until t = 108Mo. 

Figure 10 shows the scalar field as a function of areal radius and coordinate time. The scalar field begins 
radiating away soon after the formation of the apparent horizon at t = AAMq. Outside the AH, a pulse of 
scalar radiation propagates outward towards infinity. Inside the AH, the value of <j> is frozen in time because 
of the "collapse of the lapse". A large gradient in <f> develops near the matter surface just outside r s /M = 1. 
Until very late times, this region is well-covered by our grid, which follows the particles at every time step. 
However, the grid coverage in the exterior is extremely sparse: At t = 8OM0, only 18 of our 256 grid points 
lie outside coordinate radius r/Mo = 10, and at t = 90Mo, only 29 of these grid points lie outside r/Mo = 1. 
This makes it difficult to resolve the outgoing pulse of scalar radiation at late times. 

Unlike the general relativistic case discussed in section V.A, the system has not settled into its final 
state by the time the SA code loses accuracy. This is because the spacetime remains dynamical after the 
formation of a black hole — it continues to emit scalar gravitational radiation. In addition, at t = 8OM0 one 
cannot accurately determine how much energy has been lost owing to gravitational radiation because the 
tail end of the waveform has not yet propagated far enough from the black hole to justify using linearized 
theory. By t = 90Mo, the waveform has propagated a bit further into the weak-field region, but the grid 
coverage in the exterior is becoming sparse. 
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Figure 9 Oppenheimer-Snyder Collapse in Brans-Dicke theory using maximal 
slicing, and calculated by the SA method. The five solid lines represent world lines 
of matter elements containing 20%, 40%, 60%, 80%, and 100% of the interior rest 
mass. Time t and areal radius r s are measured in units of Mo, the initial mass of 
the configuration. The dotted line is the apparent horizon, which forms at about 
t = 44M . 



The SA method fails here because it cannot maintain accuracy far enough into the future. Although 
one might be able to determine the waveform and radiated energy by increasing the grid coverage by a factor 
of 4 or so, this would be extremely costly in terms of computer time. The above simulation with 256 grid 
points requires 8 hours of CPU time on a SUN Sparcstation 5 to integrate to t = 90Mo. With 512 points, 
it takes 32 CPU hours to get to t = 90Mo, and 6 additional CPU hours to reach t = lOOMo. At each time 
step, the computer spends most of its time moving the particles and calculating source terms. The length of 
each time step is small (At ~ 2 x 10 _3 Mo at late times, using 256 grid points) because of the Courant limit. 

In the AHBC method we no longer have an exponentially varying metric, so we are able to integrate 
much farther into the future. The savings in computer time are enormous: Because we no longer have a 
Courant limitation and we no longer need to move particles, the AHBC method with 256 grid points requires 
only 20 minutes of CPU time to run from t = 45Mo to t = 300Mo. The average length of each timestep is 
At ~ 0.09Mo, a factor of 45 larger than in the SA method for the same number of grid points. Furthermore, 
each step in the AHBC method takes an average of 0.4 CPU seconds, as opposed to about 1 CPU second 
for the SA method. 

In Figure 11 we show the evolution of <f> calculated using the AHBC method with 256 grid points. The 
initial data is provided by an SA slice at t = 45Mo. We show results only up to t = 216Mo, but our code 
is capable of running indefinitely. One can see the scalar field radiate away to infinity after the black hole 
forms. At the end of the simulation, we are left with the Schwarzschild solution and a constant scalar field, 
in agreement with Hawking's theorem[30] and with the simulations [4] of Shibata et ah. 

Mass conservation is shown in Figure 12. The upper two plots show M T and M s at different radii, 
calculated from Eqs. (2.64) and (2.66). Notice that values of the tensor mass measured at different radii do 
not agree, even at the initial time step. The same is true for the scalar mass, although it is difficult to see 
this from the figure because the scale of the M s plot is much larger than that of the M T plot. This small 
discrepancy is not due to numerical error; it results from second-order terms that are not accounted for in 
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Figure 10 Scalar field <f> as a function of areal radius r s on selected maximal time 
slices for the collapse shown in Figure 9. Time and distance are measured in units 
of Mo. We plot r s /Mo on a logarithmic scale in order to show the region inside the 
AH. 
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Figure 11 Same as Figure 10, except calculated using the AHBC scheme, which 
allows one to integrate farther into the future. Time and distance are measured in 
units of Mq. The areal radius r s is plotted on a linear scale. 



the linear-order expressions (2.64). These terms produce 0(1/ r) corrections to the masses. It is interesting 
to note that the instantaneous values of M T do not exhibit this error after the scalar field has radiated 
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Figure 12 Mass conservation for the AHBC method. In the upper two plots, the 
instantaneous masses from Eqs. (2.64) (solid lines) and the time-integrated masses 
from Eqs. (2.66) (dotted lines) are shown versus coordinate time at coordinate radii 
80Mo, 75Mo, 50Mo, and 25Mo. Observers at smaller radii see the outgoing pulse 
at earlier times. The lower two plots are obtained by extrapolating the data in the 
upper two plots to r = oo. Time is measured in units of Mq. 



away — the solid lines in the upper left plot converge to the same value for t > IToMq. This is because in 
the final state the spacetime is Schwarzschild, and in this case the second-order terms vanish [18]: 



A 1 ' 2 = 1 



M 
2r 



° l * 



GR, vacuum. 



(5.7) 



To remove the effects of the second-order terms and to determine the masses that one would measure 
at infinity, we extrapolate radially to r = oo by fitting to the form M = C\ + C'2/r. The result is shown 
in the lower plots of Figure 12. Because this is a radial extrapolation with no time dependence, it is only 
meaningful in the initial and final states. In the final state, the instantaneous and time-integrated values 
agree remarkably well for both M T and M s . This is a nontrivial result that verifies the accuracy of our 
numerical code. 

The final mass of the configuration is M = M T + M s ~ 0.83Mo — the system loses about 17% of its 
initial mass because of scalar gravitational radiation. While all of the scalar mass is radiated away during 
the collapse, the loss of tensor mass is small, accounting for only about 3% of the mass loss. Although the 
tensor mass should decrease monotonically with the emission of gravitational radiation[17], the instantaneous 
values in Figure 12 do not. This is another effect of the second-order terms that are not included in the 
linearized equations. 



CONCLUSION 



By adopting an AHBC method we have avoided the troublesome coordinate pathologies of singularity- 
avoiding schemes used in numerical relativity. We have used this method to evolve spherically symmetric 
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spacetimes in Brans-Dicke theory. In particular, we have treated Oppenheimer-Snyder collapse to a black 
hole in both the general relativistic limit and in a regime in which the scalar field effects are strong. We have 
also tested our code against an independent method by evolving scalar perturbations on a Schwarzschild 
background. Our code can handle black hole spacetimes that contain gravitational radiation and is capable 
of maintaining high accuracy for an arbitrarily long time. Although we have restricted ourselves to spherical 
spacetimes, our scheme does not explicitly utilize results specific to spherical symmetry (such as matching 
to an exterior Schwarzschild metric). Accordingly, we feel that generalizing the scheme to multidimensions 
holds considerable promise. 
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APPENDIX A. FINITE DIFFERENCE EQUATIONS FOR SA METHOD 

In this appendix we discuss the numerical details of the SA method, including the finite difference 
approximations for equations found in section II. The numerical grid extends from r± = to r 8max = r max . 
All variables are centered at half-grid points rj + i/2- Time steps are labeled by the index n. We define the 
time derivative operator 



T[ Y ]i + l/2 = ( Y ,t)"+l/2 

_ At„_i r i + i/2 ~~ r i + l/2 At n Y i + l/2 ~ Y i + l/2 



Y n +) n -Y ; n ,, ln At Y„." --Y" -1 



At n + At n _! At n At n + At n _! At n _! 

and a Laplacian-like operator that takes two arguments: 

( xl H0,- +1/2 

r /v n — v n \ /v n — v n N 

_ 1 x „ I »+3/2 -S- + 1/2 \ _ x „ I a + 1/2 -S-l/2 

/p3 ,p3 2 + -L I r /l r /l J * I r /l r /l 

1 i + l 'i I \'i+3/2 1 + 1/2/ \'i + l/2 ' i-\j2 , 



(Al) 



(A2) 



Here At n = t n +i — in, and X and Y are any variables defined at half-grid points rj + i/2- The quantities X" 
that appear in the finite difference expression for TZ are obtained from X" +1/2 by linear interpolation in r 2 . 

We also define an operator for a derivative with respect to r 2 : 

Q\Y]? + i,2=(Yr 3 )" +1/2 . (A3) 

This operator will only be used for variables that satisfy 

Y r 2 ~ constanti + constant2 • r 2 . (A4) 

Because all quantities are defined at half-grid points, the finite difference approximation of Q[Y]" +1/2 is 
obtained by first calculating 

( Y r2) ; ^ r r /2 " r r 1/2 (as) 

i + l/2 ' i-1/2 

at all grid points r 8 - and then interpolating linearly in r 2 to the half-grid points rj + i/2- To determine 
Q[y]" +1 12 at i = 1 we extrapolate linearly in r 2 using the values of (Y r 2j . at i = 2 and i = 3. 
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1. Wave Equation 



The finite difference form of Eqs. (3.14) can be written 



^Ktf+i/2 = 2r,- + i/2^ +1/2 QK]? + i/ 2 " < + i/ 2 n?+i/ 2 , (A6a) 

6 



T[n]? +1/2 = 2r,- +1/2 ^ +1/2 Q[n]? +1/2 - (An )3 1Z [r 3 Aa,(] n i+ll2 



(^ + i/ 2 ) 3 ( 3 + 2^)' 



(A6b) 



Notice that all quantities on the right-hand side of Eqs. (A6) are evaluated at timestep n. For a grid uniform 
in r, these equations are second-order accurate in both space and time. 

Because of the way we have finite differenced the equations, we have already included the regularity con- 
dition (3.15) at the origin, so this equation need not be imposed explicitly. A finite difference approximation 
of the boundary condition (3.17), which is imposed at the outermost grid point, is 

v n + l _ r im,i-l/2 yn , ( 1 ~ C \ I v n r »max- 1/2 v n + l \ ( « 7 \ 

{ - +1 ' 2 " ^'""- 1/J + v^Tc ) r- +1/2 " ^^ Yi -~ 1 ' 2 ) ' (A7) 

where 

C = — , (A8) 

r »max+l/2 - n max -l/2 

and Y is either II or £. This equation is second-order accurate in space and time. 

For a uniform Cartesian grid, the Courant stability condition for the above difference scheme is 

Ar 
At < t^ ;— r, (A9) 

where Ar = r;+i — r 8 -. For a variable timestep and a nonuniform radial grid, we use the stability criterion 

At = e mini \ -, —, — r \ , (A10) 

where typically we choose e = 0.5 for accuracy. 



2. Constraints 



The initial guess for <f> is computed from Eq. (3.24). We use the finite difference scheme 

= 2r,- +1/2 /3T +1/2 QM? +1/2 + /3 ' +1 r /2 ^' +1/2 " J< + i/ 2 ^ + i/2(^t)? +1/2 . (All) 

^'s + l/2 1 

Since the evolution equation is only used for an initial guess, the stability of the scheme used to solve it 
is irrelevant, and the accuracy is only important for reducing the number of iterations needed to solve the 
constraints. 

A second-order accurate finite-difference form of the momentum constraint (3.20) is 

Zj+i/2 ~ Zj-i/2 _ 8ir(S r )i _ 2 6 IF +1/2 - II 8 -_i/2 _ ( 2loiJ)1 \ 6+1/2 - 6-1/2 ( \\2) 

5 «5 £« £ ^ «2 «2 \ ^ J. I ^.2 ^.2 ' ^ ' 



r 



i + l/2 r i-l/2 5r « 5 r «' + l/2 r i-l/2 V 5< ^' / r i + l/ 2 r i-l/ 2 
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The quantities A{, IE, and <f>i defined on grid points i are determined from their values on half-grid points 
using linear interpolation in r 2 . To determine Z3/2, we integrate from r = to r = r 3 / 2 using the fact that 
t/>, £, <f>, and S r /r obey Eq. (A4) near the origin. The result is 



Z- 



V. 



3/2 



3/2' 3/2 



8^)3/2 -2r 3/2 fe 



n 



3/2 



r 5/2 r 3/2 



w n 3 / 2 ^ £5/2 - ^3/2 
^3/2 



'5/2 3/2 



(A13) 



Eqs. (A12) and (A13) are solved for Z by starting at Z3/2 and working up to ^i max +i/2- Because Eqs. (A12) 
and (A13) depend on t/> through the scalar wave terms, we must solve this equation at every step in the 
iteration of the Hamiltonian constraint. 



The Hamiltonian constraint, Eq. (3.21), can be written 

6 



^1/2 



rV/Vr= 



£)iE(*y(*), 



(A14) 



where the coefficients R^> do not depend on t/>. We solve this nonlinear equation by an iterative scheme. 
Let ip be an initial guess for t/>, and write t/> = t/>(l + (t/> — t/>)/t/>). If we have a good initial guess, we can 
evaluate the right-hand side of Eq. (A14) to first order in (t/> — VO/V 1 - The result is 



6 (rV /2 VV 2 ) 3 = ^ 1/2 J2 ^ ( * ) ^ p( * )_1 (*>(*# + (1 - P(k))i>) ■ 



(A15) 



We finite difference this equation as follows: 



6TZ 



r 3 ^/ 2 ,^ 



i+1/2 



^+1/2^2 R( i + l/2$+l/2 1 (p(k)*Pi + l/2 + (1 ~ p(k))lj>i + 1 /2 



(A16) 



This is a tridiagonal set of linear equations that can be solved for t/>. The coefficients R i+1 / 2 appearing in 
Eq. (A16) are given by 
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(A17a) 
(A17b) 
(A17c) 

(A17d) 
(A17e) 



and the values of p(k) for h = 1, . . ., 5 are —7, —1, 5, 1, and 1, respectively. 

The boundary condition at the origin, Eq. (3.22), is taken care of automatically by the finite difference 
scheme and does not need to be added explicitly. The boundary condition (3.23) is imposed at the outermost 
grid point. In finite difference form, this condition is 



(W>) 8max+1/2 - (nft),- m „_ 1/2 

r »max+l/2 - n max -l/2 



(rll) 



imax+1/2 



(rll),. 



ax-1/2 



(A18) 
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After obtaining the initial guess from Eq. (All), we iterate Eq. (A16) until convergence. Because the 
variable t/> appears in Eq. (A12), we must recalculate Z at every step in the iteration. 



3. Lapse and Shift 



Equation (3.25) is linear in a, and can be solved using the finite difference approximation 
6 



(A'+i/ 2 ) : 



:TZ[r 3 A,a} 



+ 1/2 



=«i + l/2 



£ + 1/2 \ , wII i + l/2 , 3, 



ti+i/2 \ 2 + 3/loJ <j>f +1/2 2 J 



(A19) 



This is a tridiagonal system of equations, which is easily solved by standard methods. The value of A at a 
grid point i is obtained by linear interpolation in r 2 . 

There is no need to explicitly impose a boundary condition at the origin. The boundary condition at 
infinity, Eq. (3.27), is imposed at the outermost grid point. A finite-difference version of this equation 



( ra ), m « + l/2 " ( r «X max -l/2 1 1 /, „, , ^ 



(A20) 



For the shift equation (2.24), we use the finite difference approximation 

l (h±±H - ^l!l\ = - 3a i(A' T ) 8 , 

J*i + l/2 - n-1/2 V r »' + l/2 n-1/2/ 2r; 



(A21) 



where a and K T at the grid points r 8 - are calculated from their values at the half-grid points by linear 
interpolation in r 2 . The finite difference form of the the boundary condition (3.28) is 



A m „+i/2 = 2 ( rA 'T) imax+1/2 + 2 ( rII )i max +i/2 • 
Eqs. (A21) and (A22) are easily solved by starting at i = i max and proceeding down to i = 1. 



(A22) 



APPENDIX B. INITIAL DATA FOR SA METHOD 



If one chooses Schwarzschild coordinates such that 

ds 2 = _ e 2* df 2 + e 2A dr 2 + r 2 dn2) (B1) 

then the Brans-Dicke field equations for a static configuration of dust can be written in the form[23] 

T \ l + r s t/>' /2-2e 2A \ (2 + w)r s t/>' 2 



2A' = 87rr s e 2A -^ [p + 



^'' = _^(l + e 2A) +87re 2A-^ 



3 + 2w/ 2 + r s ip' \ r 
T r s ip' 



2 + r,ip' ' 



3 + 2w 



T 



3 + 2w 



$' 



r s (2 + r s V) 



-2r s V'+-(r s V') 2 + e A -l 



(B2) 
(B3) 
(B4) 
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where 



t/> = In <f>. 



(B5) 



Here primes denote differentiation with respect to the areal radius r s , and we have chosen units such that 
<f> = 1 at r s = oo. 

For a static distribution of particles with comoving density p*, the source terms in the above equations 
are given by 



P = P*(w°) 2 , 
T = -p*. 



(B6) 
(B7) 

(B8) 

where u° is the time component of the particles' mean four- velocity field. If all particles are static, u° = 1. 
If particles are in randomly oriented circular orbits, 



Here 



u° = «°e*, 



M° = (l-r s $')" 1/2 



(B9) 



We can find the isotropic radius r from the relation 



(BIO) 



The amount of rest mass M rest enclosed within a radius r s is given by 
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Airpe r 
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(Bll) 



Near the center of the star, Eqs. (B2)-(B11) reduce to 

V = Vo + 
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(B12) 
(B13) 
(B14) 
(B15) 
(B16), 



where a zero subscript indicates a value at the center. The constants D, $o, and t/>o are determined by 
matching to the exterior solution. 

For the exterior solution, we match to the Brans type I metric[31]. This solution can be expressed 
analytically in Schwarzschild-like coordinates[32]: 



: l-£2 ■ 



„* _ cQ-x 



2( 



(B17) 
(B18) 
(B19) 
(B20) 
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The constants B, Q, and x are determined by matching to the interior solution, subject to the constraint 

Q 2 + X 2 (l+^)-xQ-l = 0. (B21) 



2 
In the exterior, the isotropic radius is given by 

1 +£ 
r = B^-^. (B22) 

Note that our parameter x corresponds to Matsuda's A/X and Brans' C/X, and our parameter Q corresponds 
to Matsuda's (1 + A)/X and Brans' (1 + C)/X. 

To determine the complete solution for a star with a given areal radius R s and a given density profile, 
we first integrate Eqs. (B2)-(B4) and (BlO)-(Bll) from r s = to r s = R s . These equations are ODEs and 
can be integrated to arbitrary accuracy by standard numerical methods. To perform the integration, we 
make an arbitrary choice for p^(r s = 0), and we set D = 1 and $ = V'o = 0. We will later determine the 
true values of D and ipo by matching to the exterior solution. The value of $o is completely arbitrary, and 
only serves to determine the time coordinate. 

After performing the integration, we match e A and if) 1 to the exterior solution, thereby obtaining ^5, 
the surface value of £, as well as the parameters Q and \. The quantity £g is given by 

ai (£| + 1) + a 2 ts = 0, (B23) 

where 

ai = e- As (l + K^'g), (B24) 

a 2 = -1 - e~ 2As (l + r s if' s + (l + |w) (r^f) • (B25) 

The constants Q and x are found from 



1 + 6 - 2j e~^ 

i-q 



q= ^? ;r . ( B26 ) 



2r s ip'sise' 



-A £ 



X = 7 S % • (B27) 

Knowledge of £g an d X determines the value of if> at the surface as calculated from the exterior solution: 

^?„ t =xln&. (B28) 

In general, this will not be equal to the value of if)g obtained from the interior solution, ipSi at i because we 
chose ip(r s = 0) = when integrating the ODE's. In order to match the value of if> at the surface, we define 



and make the transformation 



V>o = xln^-V>5 int , (B29) 



if —> t/> + ipo, 

p->pe^% (B30) 

T ->Te^">, 
M rest ^M rest e^, 
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everywhere in the interior. We transform the variables p, T, and M rest along with t/> so that the quantities A, 
$, r, r s , and ip' remain invariant — this way we do not need to recalculate £s, Q, and x from Eqs. (B23)-(B27). 
This invariance is easily verified by examining Eqs. (B2)-(B4) and (BlO)-(Bll). 

Finally, we determine the value of B from Eq. (B18) evaluated at the surface 

B=Rs (l-mt\ (B31) 

4 K J 

and we obtain the value of D by matching r s and r at the surface: 

Making the transformation r — ► Dr everywhere in the interior completes the solution. 
The various masses are found from the exterior solution: 

M s = -B X , (B33) 

M T = B(2Q - x ), (B34) 

M = 2B(Q- X ), (B35) 

Given an areal radius R s , a density profile p^(r s )/ p^(0), and a prescription for determining the four- 
velocities of particles such as Eq. (B9), the entire solution is a determined by only one parameter, the initial 
value of /?*(0) used to begin the integration (this is not the true value of p^(0) because p^(0) is modified 
by the transformation (B30)). To construct a solution with a particular mass, we vary the initial value of 
P-k(O) until Eq. (B35) yields the desired result. This procedure can be thought of as finding a root of a single 
(complicated) function of one variable. 



APPENDIX C. SOLUTION OF WAVE EQUATION IN AHBC METHOD 



In the AHBC method, the coordinate velocity of an outgoing light ray, given by Eq. (2.68), is negative 
at coordinate grid points located inside the CCH. Therefore, in this region information propagates in the 
inward direction with respect to the coordinates. In principle, any quantity at a given grid point in this 
region is completely determined by information that has propagated inward from the exterior. This implies 
that no boundary condition should be imposed at the innermost edge of the grid, since this edge lies in the 
causal future of the remainder of the spacetime. We take advantage of this property by calculating quantities 
at a given grid point inside the CCH without using information from the region interior to that grid point. 
Not only does this permit us to use a numerical scheme in which information inside the CCH cannot move 
outward, but it also frees us from imposing an explicit boundary condition at the innermost grid radius. 

Alcubierre[ll] has constructed an implicit scheme for solving the scalar wave equation in flat 1+1 
dimensional spacetime, but in an arbitrarily moving coordinate system. His scheme treates transitions 
between regions of the grid where light rays can move in only one coordinate direction, and regions in which 
they can move in both directions. In the former regions, his scheme ensures that information propagates in 
only one direction. 

We use a method based on Alcubierre's analysis, but we do not use the same finite difference scheme. 
Like his method, ours is implicit, uses spatially averaged time derivatives, and requires no rewriting of the 
field equations in an inertial coordinate system. This makes it different from the causal differencing method 
of Seidel and Suen[9] and the causal reconnection scheme of Alcubierre and Schutz[10]. However, while 
Alcubierre solves the second order wave equation for <f> using a scheme with three time levels, we solve two 
coupled first order equations (Eqs. 4.4 and 4.5) for the variables $ and n using a two-level scheme. Directly 
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solving second order equations (in time) is not very well suited for 3+1 numerical relativity, in which one 
prefers to have initial data defined on a single Cauchy surface, and one propagates this data from one time 
slice to the next. 



1. Finite Difference Equations 



We define the following time derivative operator T4 that averages over the three nearest spatial grid 
points: 



n» TT n + 1 - TT 

T [V] n + ' u ' z '- " '- ■ ' NN JJ " ; _ ' ■ " "- 1 - 1 
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(CI) 

(C2) 
(C3) 

For a spatial grid uniform in 77, we have A 8 - = 1. The above operator is second order accurate in both space 
and time. The quantity Oi is a numerical coefficient that describes the amount of spatial averaging, and 
is discussed further below. As discovered by Alcubierre, this averaging is important for solving the wave 
equation in the regime where the coordinate speed of outgoing light rays is negative. For Oi = there is no 
averaging; the time derivative at spatial grid point r\i is computed using only quantities at r\i. 

In order to average quantities in time, we introduce the operator 

Apir +i/2 = ^(^ n+i +^ n ), (C4) 

which is accurate to second order in At. Because most of our variables, such as II, are defined on the spatial 
grid points \i\, . . . , i m &x}, but the variable $ is defined on the half grid points {73/25 • • • > W + 1/2}, we 
define an operator that averages over spatial grid points: 

A [>1? + i/2=^ + i +Y")- (C5) 

Because of our staggered grid, we use a three-point spatial derivative operator T> and a two-point spatial 
derivative operator T>\j2'- 

v\y\ ee (y„). = n-*- 1 ( Y ^- Y >) + »»+i - * (Yi-Yj-A 

' * r) i+1 - 77, _i \r) i+1 - r/i J r) i+1 - 77, _i \rn - r)i-\ J 



and 



Vi f >\Y]i - tf ,).■ = ^ +1/2 "^ 1/2 - (C7) 

Vi+1/2 — Vi-1/2 



The operator T> is second order accurate even for a nonuniform spatial grid. 

In vacuum, the wave equation (4.4)-(4.5) can be written in the following implicit finite difference form 
using the above operators: 

T A [n]^ +1/2 = QV[A t [n]]- +1/2 + PVi/2[Mn- +1/2 + sA t [Av[n" +1/2 , (C8) 

+ uvvMiimnx/? + vAAA^nitlH ( c9 ) 
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The coefficients Q, P, S, W, T, U, and V are given by 

Q=^-, (ClOa) 

ri 

P=- , , ' , (ClOb) 

5 = ^( 2 + ^M? + |r^]?). ( C10c ) 

r i + l/2 

X = — !— D1/2 [C+i/2 , (ClOe) 

tr= _ ^[< +1 / 2 ) 
r i + l/2 

V=--^—T> 1/2 [aV +1/2 . (ClOg) 

Notice that the coefficients defined by Eqs. (CIO) are not centered at timestep n + 1/2, but instead are 
centered at timestep n. This is because the quantities a, /3, and t/> are only known at timestep n when the 
wave equation is solved. As a result, this difference scheme is only first order accurate in time, although it 
remains second order accurate in space (for a uniform spatial grid). However, in the case where the metric 
coefficients change much more slowly than the wave variables II and $, the scheme becomes second order 
accurate in both space and time. 

If one neglects the terms with S, X, and V in Eqs. (C8) and (C9), and one assumes a uniformly spaced 
grid, a Von Neumann analysis shows that the above difference scheme is unconditionally stable. The terms 
containing S, X, and V should not significantly affect stability because these terms do not contain derivatives 
of $ or II. 

The boundary conditions at the outer grid point, Eqs. (4.6) and (4.7), can be written 

Tjn + l Tjn i o 

1 L- 1 ^ = --vyMtinr 1 ' 2 - -A n [Mni +ll \ (en) 

LXt Ti Ti 

^Rfi + 1 /f>^ 1 9 1 

1 " '" = --Vi/Mtinr 1 ' 2 - -MMni +ll2 + -AM. (ci2) 

At ri ri ri 

Both of these conditions are imposed at i = i max . For a grid uniform in r/, they are second-order accurate in 
both space and time. 



2. Causal Solution Method 



To solve the wave equation at each time step, we first determine the location of the CCH to the nearest 
grid point using Eq. (2.68). We then define a causal boundary at r 8cB = rcB, which we place either at the 
CCH or the AH, whichever is smaller. 

Eqs. (C8) and (C9) together with the boundary conditions (Cll) and (C12) comprise a coupled system 
of linear equations for the variables $"i/ 2 an d n" + , where i = {1, . . ., i max }- Because there are a total of 
2«max variables and only 2i max — 2 equations, the system is underdetermined. However, this is a shortcoming 
of our finite difference approximation rather than a property of the underlying differential equations. In the 
continuum limit, Eqs. (4.4) and (4.5) together with the boundary conditions (4.6) and (4.7) should uniquely 
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determine $ and II everywhere in the spacetime region covered by our grid, given appropriate initial data. 
No additional boundary condition is needed at the inner edge of this region because no information can 
propagate outward from this boundary. 

This leads us to the following questions: Which of the 2i max variables should be determined by the 
2«max — 2 finite difference equations, and how should we determine the remaining variables? 



n+1 








t oooo» kf~ » ~j3 »oooo 

l n 

i=l 2 3 4 5 6 7 



n+1 



0000««»»»0000 




B. 



t n 00 

i=l 2 3 4 5 6 7 

Figure 13 Spacetime diagram showing the grid points involved in Eq. (C8) for 
i = 4. This equation is centered at the event (i = 4, t = t n +i/2)> indicated by a 
cross in the figure. The solid circles denote the grid points involved in the equation; 
other grid points are shown as open circles. In case A, left-directed and right- 
directed light rays move in opposite directions, as indicated by the light cone. In 
case B, left-directed and right-directed light rays move to the left with respect to 
the coordinates. 



For the moment, consider only Eq. (C8). For a particular value of i, this equation involves 5 grid 
points at timestep n and 5 grid points at timestep n + 1. These are the points labeled by i, i + 1, i — 1, 
i + 1/2, and i — 1/2. Figure 13A is a spacetime diagram showing these grid points for the case i = 4, in a 
coordinate system where oppositely directed photons move in opposite directions on the grid. In this case, 
one traditionally uses Eq. (C8) to determine n^" 1 " in terms of quantities defined at the other nine grid points. 

Now consider the case in Figure 13B, in which the coordinates are chosen so that both left-directed and 
right-directed light rays move to the left. We could proceed in the same way as we did in Figure 13A, and 
use Eq. (C8) for i = 4 to determine n^" 1 " . This approach should cause no difficulty for either the stability or 
accuracy of the scheme. However, we instead choose to exploit the causal structure of the problem by using 
Eq. (C8) for i = 4 to determine Ilg" 1 " . In this way, quantities at the point (i = 3,t = t n +i) only depend on 
data from points with i > 3. This would not be permitted for the case shown in Figure 13A, since in that 
case, quantities at the point (i = 3,t = t n +i) should be determined from information that propagates from 
both directions. However, in the case shown in Figure 13B, information can in principle only propagate to 
the left, a property which our scheme enforces. 

We therefore adopt the following solution method: For grid points i = 1, . . ., icB> we solve for n" + 
and $"i ; 2 using Eqs. (C8) and (C9) centered at i + 1 and i + 3/2, as in Figure 13B. For grid points 
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i = icB + 2, . . . , i max — 1> we solve for II™ + and $"i /2 using Eqs. (C8) and (C9) centered at i and i + 1/2, 
as in Figure 13A. The quantities IF max and ( I > i max +i/2 are determined by the boundary conditions (Cll) 
and (C12). Following Alcubierre, we determine the two remaining variables, IF CB _|_i and $ 8 ' CB +3/2 ) by 
requiring the functions II(r) and $(r) to be smooth: we use quadratic interpolation to obtain the two extra 
equations. 

This procedure can be encoded into a single matrix equation, which we write schematically in the form 
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^jcb + 5/2 



\ Smax + 1/2/ 



RHS. 



(C13) 



The last two rows of the matrix equation represent the boundary conditions (4.6) and (4.7), and the 
rows corresponding to i = icB + 1 and i = icB + 1/2 represent the two interpolation equations. All other 
rows represent Eq. (C8) or (C9) for some particular value of i. All quantities located on time slice t n are 
absorbed into the right-hand side (RHS) of the matrix equation. Nonzero entries in the square matrix are 
indicated by either dots or crosses. The single cross on the diagonal element of each row denotes the grid 
point being determined by that equation. The square matrix is band diagonal, so Eq. (C13) is easily solved 
by standard methods[21]. 

The above algorithm has the important property that grid points inside of the causal boundary cannot 
influence grid points in the exterior. This is true not only in the continuum limit, but in the discrete case as 
well. For example, since II^ + appears in only one equation, it must be determined by that equation; hence 
it cannot possibly affect $ or II at any other grid point. 

Until now we have not specified the spatial averaging parameter Oi that appears in the time derivative 
operator (CI). Normally one would set Oi = 0, since averaging time derivatives over space makes the 
difference scheme dispersive. However, this choice is inadequate for the grid points i = 1, . . ., icB, where we 
solve for II™ + and $"[t /2 using equations centered at i + 1 and i + 3/2. This is because the matrix elements 

multiplying II" + and $"i ; 2 ln rows i and i + 1/2 are small in magnitude compared to other elements in 
the same row, so that when one inverts the matrix to solve for these variables, one effectively sums several 
terms that nearly cancel and then divides by a small number. As a result, the matrix inversion is unstable. 
To cure this, we set Oi = 1/(2 + 2A;) inside the causal boundary r = rcB, so that there is a stronger coupling 
between neighboring spatial grid points. Outside the causal boundary, we do not need this coupling, so we 
set Oi = 0. Using different values of Oi in the exterior than in the interior causes no problem because the two 
regions of the grid are causally disconnected, even in the finite difference approximation. 
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Because we use an implicit difference scheme, there is no stability limitation on the time step. However, 
for accuracy it is useful to set the time step so that inside the causal boundary, the past directed light cone 
of a grid point at (i,t n+ i) contains the event (i + l,tn+i/2)> as in Figure 13B. In general, this requires the 
time step to be larger than the grid spacing. Notice that a Courant type condition would produce entirely 
the opposite effect: if one decreased the time step in Figure 13B sufficiently, the past light cone of the point 
[i = 3,t = t n+ i) would not contain the spacetime event (i = 4, t = t„+i/2)- From Eq. (2.68) we take the 
center of the light cone to be at Ar = —/3At, so we impose the condition 



At 



Ar 



(C14) 



at the innermost grid point. We typically choose e = 1/2. 



APPENDIX D. AHBC SOLUTION OF CONSTRAINTS 



To solve for t/> and Z, we first find the value of t/> at the AH using the evolution equation (4.12). This 
equation also provides an initial guess for t/> elsewhere. We use the finite difference approximation 



TW\? = jr (*>M? + \w 
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where the operator T> is defined in Eq. (C6) and the operator T is given by Eq. (Al). This scheme is 
second-order accurate in space and time, even for a nonuniform grid or for unequal timesteps. The stability 
of the scheme is irrelevant since the result is only retained at the AH — at all other grid points, t/> is refined 
using the Hamiltonian constraint. 

Next, we solve the momentum and Hamiltonian constraints simultaneously using an iterative scheme. 
These equations are coupled because Brans-Dicke scalar radiation terms containing t/> appear in the momen- 
tum constraint (4.8). Let ip be an initial guess for t/>, and let Z be an initial guess for Z. If we substitute 
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into the constraint equations (4.8) and (4.9), and expand to first order in the small quantities (t/> — ij))/ij) 
and (Z — Z)/Z, the result is 
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Applying the same linearization procedure to the boundary condition at the AH, Eq. (4.10), we obtain 
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at AH. 



(D5) 



For the two other boundary conditions, we set the t/> at the AH to the value obtained from the evolution 
equation, and we use Eq. (4.11) at the outer grid point. 

The constraints, together with the boundary conditions that they must satisfy, can be written in the 
finite difference form 
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Here the operator T> 2 is defined by 
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and the operators T>, 2>i/2> an d A v are defined as in Eqs. (C5)-(C7). In the above equations, the averaging 
operator A v takes precedence over other operations, e. g., 



A V [XY 2 ] = A V [X](A V [Y}) 2 . 



(D12) 



Note that £>[£>[Y]] # T> 2 [Y], but for an equally spaced grid, 



Vi,2[D 1 ,2[Y]}=V 2 [Y}. 
43 



(D13) 



Eqs. (D6)-(D10) can be encoded into a single matrix equation: 



x • • • 
x • • 



x • • • 

x • • • 

• x • • 

x 

• • x • 

• • x 



\ / \ 



• • x • 



V 



x/ 



^AH-1 


V'iAH-l 


'"•AH 


V^AH 


^•'ah + 1 


V'iAH + l 


2, m „ 


V </>*_ ) 



RHS. 



(D14) 



Here crosses denote elements on the diagonal, and dots represent all other nonzero entries. The RHS matrix 
does not depend on t/> or Z , but does depend on the initial guesses ip and Z. 

Given an initial guess, we iterate the matrix equation (D14) until convergence, using the values of t/> and 
Z at each step as initial guesses for the next step. The matrix equation is solved by a standard band-diagonal 
inversion technique [21]. 



APPENDIX E. AHBC SOLUTION OF LAPSE AND SHIFT EQUATIONS 



The lapse and shift equations, Eqs. (4.13) and (4.14), together with the boundary conditions (4.16)- 
(4.18), can be written in the following finite difference form: 



V 2 [a]i + V[a] 



1 + 2D[i/>]j 



? (( * T)lw *)' + $ *+* 



fi + vMiW] 


V i 
z 

r 3 t/> 6 


i J J . 

= 
i + 1/2 



!^ + £(2> 1/2 [*].- + 2*,. 



'E > i/2[X]i+i/2 + y-A v [a]i+i/2-A v 



X>i/ 2 [ar] i _ 1 /2 = 1 + - (Uin + n 8 _ir 8 _i) , 

r i-l/2 & 



Xi 



(K T )i iii 



2 ' 



(El) 


(E2) 


(E3) 


(E4) 


(E5) 
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Here X{ = fii/ri, and 



(^i).- 



(F2), 



yb r ji 



Mi 



t. 



2 / ' 



8tt ns r r )i Tj (5 r ),- 

^ # 3 + 2w + V? 



<«X^<n*f-*>--^ + ^-*-*™ 



nfaipi 



(jr 4) ,.= (^? W n,. _$,.)• 



(E6a) 

(E6b) 

(E6c) 
(E6d) 



The operators ^, £>i/2, T>, and £> 2 are defined by Eqs. (C5)-(C7) and (Dll). The averaging operator A v 
has precedence over other operations. 

The lapse and shift equations are coupled because of the boundary condition (E3). To solve them 
simultaneously, we write Eqs. (E1)-(E5) as a single matrix equation: 



/ 



x • • • 
x • 



\ / \ 

Xr 



x • • • 



X • • • 
• X 

X • • • 

• X • 



V 



x/ 



• ah- 



a; 



AH — 1 

Xi 



UH 



a. 



AH 



^AH + 1 
a >AH + l 

\ «Wx / 



RHS. 



(E7) 



Here crosses denote elements on the diagonal, and dots represent all other nonzero elements. This equation 
is solved for a and X by a standard band-diagonal inversion method [21]. 
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